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Abstract. In this paper we study pattern formation arising in a system of a 
single reaction-diffusion equation coupled with subsystem of ordinary differen- 
tial equations, describing spatially-distributed growth of clonal populations of 
precancerous cells, whose proliferation is controled by growth factors diffusing 
in the extracellular medium and binding to the cell surface. We extend the 
results on the existence of nonhomogenous stationary solutions obtained in [9] 
to a general Hill- type production function and full parameter set. Using spec- 
tral analysis and perturbation theory we derive conditions for the linearized 
stability of such spatial patterns. 

1. Introduction 

Partial differential equations of diffusion type have long served to model regula- 
tory feedbacks and pattern formation in aggregates of living cells. Classical mathe- 
matical models of pattern formation in cell assemblies have been constructed using 
reaction-diffusion equations. They have been applied to describe pattern formation 
of animal coat markings, bacterial and cellular growth patterns, tumor growth and 
tissue development, see e.g., [10] and [11] and references therein. One of the mecha- 
nisms of pattern formation in reaction-diffusion systems, prevalent in the modeling 
literature since the seminal paper of Allan Turing [14], is diffusion driven instability 
(Turing- type instability). 

Diffusion-driven instability arises in a reaction-diffusion system, when there ex- 
ists a spatially homogeneous solution, which is asymptotically stable in the sense 
of linearized stability in the space of constant functions, but it is unstable with re- 
spect to spatially inhomogeneous perturbation. The majority of theoretical studies 
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in theory of pattern formation focus on the analysis of the systems of two or more 
reaction-diffusion equations. In many biological applications it is relevant to con- 
sider systems consisting of a single reaction-diffusion equation coupled with a system 
of ordinary differential equations. Such models can also exhibit diffusion-driven in- 
stability. However, they are very different from classical Turing-type models and the 
spatial structure of the pattern emerging from the destabilisation of the spatially ho- 
mogeneous steady state cannot be concluded from a linear stability analysis [8, 10]. 
The models exhibit qualitatively new patterns of behavior of solutions, including, 
in some cases, a strong dependence of the emerging pattern on initial conditions 
and quasi-stability followed by rapid growth of solutions [8] . 

Mathematical theory exists only for special cases of such solutions arising in 
the models, which can be simplified to the single reaction-diffusion equation with 
nonlocal terms and with fast growing kinetics. The example of such systems are 
Gierer-Meinhard and Gray-Scott models. The existence and stability of the solu- 
tions of such systems were intensively studied using singular perturbation analysis 
and spectral analysis of the eigenvalue problem associated with the linearization 
around the "spike- like" solution, e.g. [3, 4, 15]. The structure of the kinetics in- 
volved in the cell proliferation model considered by us is different, in particular 
the autocatalysis is excluded by the biological assumptions and the solutions of 
our model are uniformly bounded. As shown in [7, 9] the models have stationary 
solutions of periodic type, the maxima and minima of which may be of the spike or 
plateau type. Numerical simulations show that in some cases solutions of the model 
converge to a spatially heterogeneous pattern, which persists for an arbitrary long 
time while in other cases the transient growth of nonconstant pattern is observed 
and ultimately the solution converges to a stable spatially homogeneous state. 

In this paper we approach the issue of the stability of nonconstant stationary 
solutions and investigate linear stability of spatially heterogeneous solutions of the 
model of the reaction-diffusion equation coupled with two ordinary differential equa- 
tions, which was proposed in [9]. The paper is organized as follow: In Section 2 
the model is introduced and the results on existence, regularity and boundedness of 
the model solutions are presented. In Section 3 existence of a spatially nonconstant 
steady state is shown. Section 4 is devoted to the linearized stability analysis. 

2. Model description 
We consider a model of a cell population controlled by a diffusive growth factor, 

dtc = — — — d c c + (i in (0, oo) x (0, 1), 

6 + c 



(1) 



d t b = ac g — d^b — vb in (0, oo) x (0, 1), 

1 (3c k 
dtg = -d%g- ac 2 g - d g g + + vb in (0, oo) x (0, 1), 

7 1 + c R 

d x g(t, 0) = d x g(t, 1) - in(0,oo), 



with initial conditions 

c(0,x) = c in (x), 6(0, x) = b in (x), g(Q,x) = g in (x) for x 6(0,1), 

where c denotes the concentration of precancerous cells, whose proliferation rate is 
reduced by cell crowding but enhanced in a paracrine manner by a hypothetical 
biomolecular growth factor 6 bound to cells. Free growth factor g is secreted by 
the cells, then it diffuses among cells with diffusion constant 1/7, and binds to cell 
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membrane receptors becoming the bound factor b. Then, it dissociates at a rate v, 
returning to the 6-pool. Parameter fj, denotes a small influx of new precancerous 
cells due to mutation. All coefficients in the system (f ) are assumed to be constant 
and positive, and k £ N. For k larger than f , production of growth factor molecules 
is given by a Hill function and models a process with fast the saturation effect. 

Theorem 2.1. For nonnegative initial data (ci n ,bi n , gi n ) £ C 2+a ([0, f]) 3 there ex- 
ists a nonnegative global solution of system (1), (c,b,g) £ C 1+Q / 2:2+Q ([0, oo) x 

[o,i]) 3 . 

Proof. Using the existence and regularity theory for systems of parabolic and ordi- 
nary differential equations (see [5], [12]), for (cf n , b in ,g in ) £ C 2+a ([0, l]) 3 , we obtain, 
due to local Lipschitz continuity of reaction terms in the system, the existence of a 
local solution of (1), (c,b,g) £ C* 1+a / 2 - 2+Q ([0, T ] x [0,1]) 3 for some T < oo. 

The theory of bounded invariant rectangles (see [2], [13]) and the properties 
of functions F c — ^ — d c c + /i, Ff, = ac 2 g — d^b — i/b, F g = j^jr — ac 2 g — 
d g g + vb, i.e. are smooth for c > 0, b > 0, g > and in {(c,b,g) £ M 3 : c ^ 
-6}, F c (£,b,g) > for all b > 0, g > 0, F b (c,0,g) > for all c > g > 
0, and F g (c,b,0) > for all c > b > 0, imply that the set S = {(c,b,g) : 
c(t,x) > mm{mxn x& [ ^c in (x),-^},b(t,x) > 0,g(t,x) > 0,x £ [0,1], t > 0} is 
positive invariant and solutions are nonnegative for nonnegative initial data. Then, 
from the first equation in (1) follows the estimate 

d t c < (9 — d c )c + fi, 

and, by Gronwall inequality, we obtain, that c is bounded for all finite time, i.e., 
sup c(t, x) < ^exp ((9 — d c )t). The boundedness of c implies also the boundedness 

see [0,1] 

of b and g, and existence of a global solution. □ 

3. Existence of positive non-constant steady states 
Stationary solutions of (1) can be found from the system 



(2) 



9bc 

— d c c + /i = 0, 

o + c 

ac 2 <7 — df,b — i/6 = 0, 

- g" - ac 2 g - d g g + + ^6 = 0, 5 '(0) = g'(l) = 0. 

7 i + c 



We look for a positive solution (c(x), b{x),g(x)) defined in (0, 1) that has at least one 
non-constant component. It is clear that for any solution of (2) either all functions 
c, b and g are constant or all depend on x effectively. The next observation is that 
the function c doesn't change the sign. In fact, if c{xq) = for some Xq £ (0, 1), 
then it follows from the first equation that fi = 0, which is impossible. Let us solve 
the system 

( 9b c 

(3) \bTc- dcC + ^ = ^ 

I ac 2 g — d\jb — vb = 

with respect to c and b. 

Lemma 3.1. If 9 < d c , then there exists a unique solution (c(g),b(g)) of (3), 
which is continuous and positive for all g > 0. In the case 9 — d c there is only one 
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positive solution and for 9 > d c there are two positive solutions that exist only in 
some interval (0,g*), where g* is a finite number. 

Proof. Let us introduce the temporary notation £(g) = J^+ v g- First observe that 

(4) b(g) = ^£(g)c 2 (g), 

and hence b is positive as soon as g is positive. The main question is whether there 
in fact exists a positive c. Substituting expression (4) into the first equation of (3) 
yields 

(5) {9 - d c )£c 2 + fi(£- d c )c + fi 2 = 0. 

This quadratic equation with respect to c admits real roots provided 

D = fj, 2 {{I - d c ) 2 - 4(9 - d c )£) = v 2 ((£ + d c ) 2 - 4 j d c ?j > 0. 

Note that (p + q) 2 — 4epq > for all positive p, q if and only if e < 1. Hence, under 
the assumption 6 < d c , equation (5) has the solution 



£{g) ~d c + ^/(£(g)-d c ) 2 +4(d c -9)£(g) 

(6) C(5) = ^ 2 { d c -9)£ { g) ' 

which is positive for all g > 0. Furthermore, 

(7) c(0+) = £' Ci+00)= d^9- 

In the the critical case 9 = d c there exists a unique solution of (5) 

(8) c(g) = - 11 



where c is positive for g £ 0, ^^"^ J ■ In the case 9 > d c there are two positive 
solutions of (5), defined for g G (0,3*) with g* = ^-(VO - ^9-d c ) 2 , g* < 

^ctii • on ly one °f them is bounded at g = 0, and this solution is given by 
(6). Moreover the solution is also bounded at the point g*: 



(9) c(g*) 



2{Q-d c )l{g*y 



Note that g* is the smallest positive point in which the discriminant D = D(g) 
vanishes. Observe also that £(g) < d c for all g < g* . □ 

We next come back to system (2). Let u) — . In view of Lemma 3.1 the 
function g must be a positive solution of the nonlinear boundary value problem 

do + ^f>)-:n% il,(0 ' 1) ' 

U(o) = o, 9 '(i) = o, 

where c is given by (6), if 9 ^ d c , and by (8) otherwise. We must keep in mind that 
the function c depends on parameters a, /i, v 1 9, and d c . 



STABILITY OF NONCONSTANT STATIONARY SOLUTIONS 



5 



3.1. Nonlinear boundary value problem. Let us consider the boundary value 
problem 

(11) 9" = lh(g), xe(0,l), </(0) = 0, g'(l) = 0, 

where 

2 /^\ /?c fc (5) "4 



%) -rf 9 .g(i + ^c 2 ( 5 )) 



i + c k { g y dg{d b + vY 

and c is given by (6) or (8). The equation of this type, so called system with 
one degree of freedom, was intensively studied in classical mechanics. For a deeper 
analysis of the system we refer to [1, Sec. 12]. 




Figure 1. Plots of the energy U and the trajectories of a dyna- 
mical system. 



Theorem 3.2. If function h — h(g) has no less than three positive roots, then there 
exists a set T C K+ of diffusion constants 7 for which boundary value problem (11) 
admits a positive solution. 

Proof. We assume for a while that 7 = 1. Suppose that the potential energy 

U(g) = - [ 9 h(0dC 

Jo 

has a positive critical point go that is a local minimum of U. The trajectories of 
dynamic system g' = z, z' = h(g) resemble ellipses near the point (go , 0) of the 
phase space. Let us consider the trajectory £ that intersects the axis z = at 
positive points qi and 92 (see Fig. 1). Therefore there exists a periodic solution 
g = g(x), z = z(x) of the dynamic system subject to initial conditions g(0) = q±, 
z(0) = 0, and its period is given by 

dr/ 

^2(u( qi )-wwy 
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Remark that the energy U without local minimum points is either monotonic or a 
function with a single maximum point. In both cases no trajectory starting at a 
point (51, 0) can return again to the line z = 0, with the exception of the equilibrium 
positions. 

Next, g is a positive solution of the initial problem g" = h(g), g(0) = qi, g'(0) = 0. 
Moreover, g'i 11 ^-) = for any natural n, because g'i 11 ^-) = z(0) — for even n and 
g'(^) = z (^) — for odd n. Given n, we consider the function u{x) = g{^j=) that 

is obviously the solution to equation u" = jh(u), and u'(0) = 0. We set j n = n2T2 , 
so that 

1 ,/ 1 \ 1 , [nT s 



u'(l) = — g> — )=-^g' =0 - 

Hence, u is a positive solution to (11). 

It follows that any close trajectory £ lying in the half-plane g > produces 
a countable set of positive solutions u n ((, •) to (11) with 7 = 7n(C)- All these 
sequences {7n(C)}«eN form the set T, which is in general uncountable. 

The first and primary question is whether the energy U has a local minimum. 
We consider first the case 9 < d c . The energy has a local minimum at point g = go 
if only h(go) — 0, h > in a left-side neighborhood of go and h < in a right-side 
one. Let 

/ ll ( 9 ) = rf 9ff (l + u;c 2 ( ff )), ft 2 ( g ) = i ^L . 

We look for a root 50 of the equation h\{g) = h 2 (g) for which /ii(<?) > h 2 (g) the 
left of go and < ^2(9) on its right. A trivial verification shows that c = c(g) 

(given by (6)) is a steadily increasing function and 

(12) c ( +0 ) = £, c (g)^-t— as 5 ^+^. 

Set c m i n = and c max = d . It is convenient to consider the c-representation of 
functions hj on interval (c m i n ,c max ): 

, , N a(c- c min )(l +uc 2 ) /3c fc 
/11(c) = ? r , h 2 (c) - 



c(c max - c) ' 1 + c k ' 

where a = ^M~gY^ ■ Both functions are strictly increasing in the region under 

study. In regard to hi, its derivative can be written as 

, , 1 v (1 + wc 2 )(c min (c max - c) + c(c - Cmin)) + 2cjc 2 (c - c min )(c max - c) 
MC) = c^-cf 

and is obviously positive for c € (c m i n ,c max ). The function hi has two vertical 

asymptotes c — and c = c max , while /12 is bounded. In addition, h 2 has an 

inflection point c = y f^pj, whereas the inflection point of hi depends on parameters 
of the model. Fig. 2 shows the typical plots of hi and h 2 . Evidently, the energy U 
has a local minimum if plots of the functions intersect at three points. 
Similarly, for 9 — d c we obtain that c{g) is monoton increasing and 

jj, . . d c (d b + v) 



(+0) = - , c{g) -4-00 as g -> 



c> . , - 

The function /ii is nearly linear for small g and growth quadratically as g 

whereas /i 2 ( ff ) ~ cg k for 5 € (0, d ^^ 1 ) such that c( 5 ) < {/|^, and is 

bounded for g dc{d *+ v) ; additionally, hi(0) = < h 2 (0) - — — - 
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k > 2 there exist sets of parameters d c , db, v, a, \x, (3, such that the function h(g) has 

d c (d b +u) \ 



three roots in (0, 





Figure 2. Sample plots of hi and h 2 corresponding to nonexis- 
tence and existence of a nonhomogeneous stationary solution in 
the left and right panels, respectively. 



For > d c and c(g) defined by (6), we obtain 



:(+0) = 



c(g*) 



n(d c -l(g*)) 
2(e-d e )l(g*) 



and again hi is nearly linear, /12(g) ~ cg k for small g, hi(0) = < /^(O), and both 
functions are bounded on (0,g*). For k > 2 there exist two roots, as functions of 
parameters, of h(g) on the interval (0, g*). The third root of h(g) lies on the interval 
(g**,oo), where g** is the largest positive point in which D = D(g) vanishes. □ 

Remark 1. Let d c < 9, the solution c(g) of (5) is defined by 



(13) 



c{g) = m ■ 



£(g) -dc- ^(£(g)-d c ) 2 +A(d c -9)£(g) 



2(d c - 6)£{g) 



parameters d c , db, (Ji, v, a satisfy the assumption < g* < (db + v)d c /(afi), and 
(3 is choosen such that h(g*) < 0. Then the function h is continuous in (0,g*], 
and h(g) — > +00 as g — > +0. This properties of h and the assumption h(g*) < 
provide existence of a local minimum of the energy U at a point go G (0, g*). Thus, 
there exists a set r C M+ of diffusion constants 7 for which the boundary value 
problem (11), with c(g) defined by (13), admits a positive solution. This situation 
was considered in [9] with k = 1. 

If U has a local minimum go, then there are a local maximum g~ in the interval 
(0,go) and a local maximum g + in the interval (go, 00) as shown in Fig. 1. Set 
£0 = min{|.g - g~\, \g a - 9 + \}- 

Corollary 1. For any e £ (0,£o) there exists a countable set of solutions g n to 
boundary value problem (11) with 7 = j n and j n — > as n — > 00, such that 

\9n(x) - go I < £ 

for all x G (0, 1). Moreover g n is rapidly oscillating function as j n goes to 0. 
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FIGURE 3. Spatially nonhomogeneous solutions of the stationary 
problem for the different values of 7. Parameters: a = 0.1, (3 = 10, 
v = 0.1, (i = 0.1, d b = d g = 0.1, d c = 0.6, = 0.59, and 7 = 
12188.16 and 7 = 4900 in the left and right panels, respectively. 



4. Stability and destabilisation of steady states 

Let X be a complex Banach space with norm || ■ || and let A be a closed operator 
on X with a dense domain T>(A). Suppose that —A is a sectorial operator and 
Re A < for all A G cr(A). Hereafter, a(T) stands for the spectrum of an operator 
T. For each s > Owe introduce the interpolation space X s — T>((—A) s ) equipped 
with the norm ||x|| s = || (— A) s x|| . Clearly A" = X. We consider the equation 

(14) d t u = Au + f{u) 

in the Banach space X. Let u* G D(A) be a stationary solution of (14) . 

To study stability of the stationary solutions we apply the following proposition 
about stability and instability by the linear approximation that is a version of 
Theorems 5.1.1 and 5.1.3 in [5] adapted for our purposes. 

Proposition 1. Let f : hi — > X be a locally Lipschitz continuous map in a neigh- 
borhood hi C X s of a steady state u* , for some s G (0, 1). Suppose that for z G X s 
the map f admits the representation 

(15) f(u*+z) = f(u*)+Bz+p(u*,z), 

provided \\z\\ s is small enough. If 

o B is a linear bounded operator from X s to X , 
o ||p(u*,*)||=o(||z|| a ) a S ||z|| s ^0 ; 

o the spectrum of A + B lies in the set {A G C : Re A < h} for some h < 0, 

then the equilibrium solution u* of (14) is asymptotically stable in X s . Moreover, 
if the spectrum o~(A + B) and the half-plane {A G C: Re A > 0} have a non-empty 
intersection, then u* is unstable. 

4.1. Linearised analysis. For simplicity of notation, we denote the vector (c, b, g) 
by u = (ui, u 2 , U3) so that system (1) reads 



(16) 



d t u = Au + f(u) 
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in the Hilbert space X = L 2 (0, 1) © L 2 {0, 1) © L 2 (0, 1), where 

/ 9u\u 2 

' 

-(db + v) 
iV 



A = 



u\ + u 2 
au\u 3 



j- + z^u 2 — au\U3 



) 



(IgV 



and iVo is the Sturm-Liouville operator given by 

1 d 2 v 

7 dec" 1 

on the interval [0, 1], subject to the Neumann boundary conditions, 

V(N )={v€ W 2 2 (0,1): t/(0) = 0y(l)=0}. 

Recall that any self-adjoint and bounded from below operator is sectorial, [5, 
p. 19]. The operator A in X with the domain V{A) = L 2 {0, 1) © i 2 (0, 1) © V(N ) 
is self-adjoint. Moreover, —A is bounded from below, since 

a(A) = {-d c , -d b -u}U {-d g - 7~ Vj 2 }^ - 

Hence —A is a sectorial operator and Re A < for all A £ &(A). 

We can introduce interpolation spaces X s — T)((—A) s ), each of which is a Hilbert 
subspaceofL 2 (0,l)fflL 2 (0,l)©W^(0,l). Set \\p\\ s = || Pl || La + ||p 2 || i2 + ||p 3 ||wf- 

The function / is smooth in R\ = {y £ R 3 : y k > 0, k = 1,2,3}. Therefore, / 
admits the representation 

f{y + z) = f(y) + B(y)z+p(y,z), 

where the remainder satisfies the estimate 

(17) \\p( y ,z)\\ R3 <$(y)\\z\\ 2 R3 

in a neighborhood of any point y £ M 3 , with a continuous function Here 



B(y) 



(j/i + 2/2) : 

2ayi2/3 
fe-i 



{yi + yiY 
o 



o 



\ 



k \ 2 



uyl 

-ay 2 



Suppose now that u* is a positive steady state, that is Au* + f(u*) — and 
u*(x) £ R\ for all x £ (0,1). Assume also that u* are C 2 -functions, j = 1,2,3. 
Obviously, B(u*) is a linear bounded operator from X s to X, for each s £ (0, 1). 

Next, using (17) we obtain f(u* + z) = f(u*)+B(u*)z+p(u* , z) with the estimate 
of the remainder 



i|p(u*,z)|| R 3 < p{u*(x)) \\z(t,x)\\^ 3 < max p(u* (x)) \\z(t,x) \ 



< 



xe[o,i] 



Ml 



Consequently, for any z £ X s we conclude 



||p( U *,z)||< Cl (||z 1 ||| 2 + ||z 2 ||^ 



MIL 



< dOMi, + \\z 2 \\l 2 + Wz.W^) < c 2 \\z\\l = o(\\z\\ s ), 
as H^lls — > 0, with constants c k being independent on t. 
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For nonlinear model (1) the linearized system at the steady state u* can be 
written as dtz = (A + B{u*))z. Then turning back to the previous notation, we 
obtain 

obi . \ oa 



(18) 



-d, 



d t z 2 — 2ac*g* z\ - {d b + v) z 2 + acl z 3 , 



Z2, 



d t z 3 = 



(1 



„k\2 



2ac*g* I z\ + v z 2 



7 



Z3- 



The first and primary question is whether the spectrum of A + B(u*) lies in the 
half-plane {A € C : Re A < h} for some h < 0. 

4.2. Spectral analysis of linearised problem. Let us consider the Sturm-Liouville 
operator N = ±£2 - a<%{x) - d g in the domain V{N) = {v € W^O, 1): v'(0) 
0. v'(l) = 0}. For abbreviation, we write A instead of A + B(u*), and consider the 
operator 

l(x) m{x) 
n{x) —yc 
,r(x) v 



A 



with 



I = 



Obi 



2ac ir g* 

)„2 



(b. 



(1 + cJ) 2 



Of course, constant x and functions to, n, q are positive in [0, 1], since the stationary 
solution (c* , b* , g* ) is positive. 

The matrix A is regarded as an unbounded operator in the space X with the 
domain V(A) = L 2 (0, 1) © L 2 (0, 1) ffi V(N). It can be written in the form A = 
An + A\ , where 



An = 




and A\ = 




Notice that A\ is a bounded operator in X. The operator An is the direct sum 
M © N, where the matrix operator 

' I TO 
,0 -X, 



A/ = 



is bounded in Y = L 2 (0, 1)©L 2 (0, 1). Moreover An is a selfadjoint operator. Clearly 
a(A Q ) = a{M) U a(N) C R. 

In order to describe the spectrum of M we consider the equation {M — \E)cf) = f 
that is to say 



(19) 



(I - A) 01 + TO02 = /l, 
-(x + A)0 2 = /2 



for each / € Y. If A = — x, then system (19) is in general unsolvable. Hence, 
— x € a(M). In the case A ^ — x the solution of (19) can be written in the explicit 
form 



(20) 



4>i{x, A) 



m(x)f 2 {x) 



l(x)-X (x + X)(l(x) - A)' 



02 (x, A) = - 



x + A' 
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Nevertheless, this solution does not always belong to Y. Indeed, if 1(xq) = Ao for 
some point (x , Ao) G [0, 1] x K, then the function (l(x) — Ao) -1 is not quadratic 
integrable in the neighborhood of Xq, because l(x) — Ao ~ k(x — x ) s as x — > Xq and 
6 > 1 due to the smoothness of I. If l(x) — Ao is different from zero for all x G [0, 1], 
then (20) is a unique L 2 -solution to (19) for all / G Y, and so Ao is a point of the 
resolvent set of M. 

Proposition 2. The spectrum of Ao is real and consists of the point — x, the range 
of function I and all eigenvalues of N . 

Theorem 4.1. Let ft* = max {— x, Ai, max^^i] K x )}> where Ai is the largest 
eigenvalue of N . The spectrum of Aq lies in the set {A G C: Re A < h*} and 
h* < 0. 

Proof. The spectrum of Ao coincides with the set 

{-x} U {Ao}£Li U [ min l(x), max l(x)], 
J xe[o,i] see [o,i] 

where are eigenvalues of the Sturm-Liouville operator N = — acl(x) — d g 
subject to the Neumann boundary conditions. Since etc 2 + d g > in [0,1], all 
eigenvalues of N are strictly negative, and hence a(N) C {A G C: Re A < Ai}, 
where Ai < 0. Moreover, the first equation in (2) yields 

r<m u \ - 9b *^ h ^ 9K{x) rf - <n 

[ ' {X> ~ (b*(x) + c*(x)y c< h(x) + c*(x) c ~ c*(x) <U - 

Hence h* is a negative number. □ 

Next, we will consider the operator A as the perturbation of Ao by the bounded 
operator A\. Let T and S be operators with the same domain space TL such that 
V{T) C V{S) and 

(22) ||5 U ||<a|H| + 6||Tu||, u G T>(T), 

where a, b are nonnegative constants. Then, we say that S is relatively bounded 
with respect to T or simply T-bounded. Assume that T is closed and there exists 
a bounded operator T^ 1 , and S is T-bounded with constants a, b satisfying the 
inequality 

aWT- 1 ]] +b < 1. 

Then, T + S is a closed and bounded invertible operator [G, Th.1.16]. 

Theorem 4.2. Under the assumptions of Theorem 4-1 the spectrum of A lies in 

the set {A G C : Re A < h} with h < provided 

(23) 

f c^ 1 (x) I 

max^ v, a max cJx), 4a max (cJx)g*(x)),k8 max -— — — , -. NN „ > < \h*\. 
\ xe[o,i] x€[o.i] y v Jy v " *e[o,i] {l + c k ;{x)) 2 J 1 1 

Proof. Let ft, G (/i*,0). Choose A G {A G C: Re A > ft,}. Due to Theorem 4.1 the 
operator _4o — Ao-E is bounded invertible. On the other hand, A\ is bounded in X 
and 

\\AmW < ||(|r|+n)0i||£2 + ||^ 2 || L 2 + ||# 3 |U= < a\\u\\, 

where 

max (|r(x)| +n(x)), v, max \q(x)\ 
xe[o,i] x£[o,i] 
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Thus Ai is (Ao — Ao-E)-bounded with b — in the corresponding inequality (22). 
Consequently if a||(.4o — Aoi?) — 1 1 < 1 ; then the operator A—X^E = An+A\ — X a E 
is bounded invertible. Note that 

\\(Ao X E)-'\\ = dist{Xo ] a{Ao)) < j^—f^, 

since Ao is self-adjoint. If (23) holds, then a < \h*\ since |r(a;)| + n{x) = r + n = 

(i+c*) 2 ^ or r ( a: ) — ^ an< ^ l r ( x )l + n ( x ) — ^ n ( x ) f° r r i x ) < 0, x € [0, 1], q — ac\ and 
v <]h*\. Therefore a < \h* - h\ for some h € {h*,0) and a\\(A - A ^) — 1 1| < 1- 
We showed that the set {A e C : Re A > h} belongs to the resolvent set of A, which 
completes the proof. □ 

4.3. Linearized stability conditions. Theorem 4.2 implies that under the as- 
sumptions (23) the spectrum of linearized operator A lies in {A S C : ReA < h} for 
some h < 0. Then, from Proposition 1 follows the linearized stability of a stationary 
solution (c*,6*,2*) in X s C i 2 (0, 1) ® L 2 (0, 1) ® iy| s (0, 1), s€ (0,1). 

Corollary 2. A spatially nonconstant stationary solution (c*,6*,g*) of system (1) 
is linearly stable in X s if the following conditions are fulfilled 

ml(x)<\h*\, W <\h% 4a{c4x)g*(x))<\h*\, v<\h*\ 

for all x € [0, 1], where h* is defined in Theorem 1^.1. 

Remark 2. It is enough to verify (24) in a local minimum go of the potential U 
from Theorem 3.2 and for the sets of parameters such that the local minimum exists. 
Then, Corollary 1 and the strict inequality in (24) provide stability of a stationary 
solution with the small amplitude. 



5. Conclusions 

In this paper we performed linearized stability analysis of the nonhomogeneous 
stationary solutions of a single reaction-diffusion equation coupled to ordinary dif- 
ferential equations. We focused on a specific model of pattern formation in a system 
of cells with the proliferation regulated by a diffusive growth factor. We extended 
results on the stationary problem obtained in [9] to the larger space of parameters 
and showed existence of inifinitely many stationary solutions. Numerical calcula- 
tions indicate that the relationship between d c and 8 is significant for the existence 
of nonconstant stationary solutions. It is observed that for d c ^ 9 the system (1) 
admits the nonconstant stationary solutions when the difference \d c — 9\ is small 
enough. For d c — 9 the value of the parameter /3 plays an important role in existence 
of the stationary solutions. 

Performing spectral analysis using sectorial operators and the perturbation theo- 
ry we obtained conditions for the linearized stability of the spatial patterns. For- 
mulation of the conditions in terms of the model parameters independent of the 
values of stationary solutions is difficult to obtain due to the complex structure of 
the stationary problem. For some parameter values, due to existence of multiple 
solutions of the ODEs subsystem, the model exhibits also discontinuous stationary 
solutions and it is difficult to decide which patterns we see in numerical simulations. 
The mechanism of pattern selection is a topic of further analytical and numerical 
studies. 
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